; FFT program by Piotr Nowak
; 10/21/2009
; Based on FFT C code by Tom Roberts, found on               ;http://www.jjj.de/fft/fftpage.html



; fr[] and fi[] are real and imaginary arrays for input and output
; size of each array is 2^m

; mem(1000) = m
; mem(1001) = fr[0]
; mem(1001+2^m) = fi[0]
; mem(1001+2*2^m) = sinewave

bez r0, start

FIX_MPY: ; Fixed-point multiplication subroutine
; parameters passed in r12 and r13, returned in r12
; r13, r14 and r15 will be used for temporaties
; r6 contains return address
mv  r14, r12
mul r14, r13 ; r14 = lo, r15 = hi

mv  r12, r15
sll r12, 1
srl r14, 30
li r15, 1

mv  r13, r15
and r13, r14
add r12, r13
srl r14, 1
add r12, r14

jalr r0, r6 ; return

start:
li r1, 1000 ; 1000, base mem address

; Find n = 2^m
load r2, 0(r1) ; r2 = m
addi r3, 1 ; r3 = 1
sllv r3, r2

 ; nloop:
; bez r4, nloopdone
; sll r3, 1
; addi r4, -1
; bez r0, nloop;
; nloopdone: ; now r3 = n

li r4, 0 ; r4 = mr = 0
mv  r5, r3
addi r5, -1 ; r5 = nn = n - 1

; do decimation-in-time


li r2, 1 ; r2 = m = 1
ditForLoop:
mv  r14, r5
sub r14, r2 ; r14 = nn - m
blz r14, endDitForLoop
mv  r6, r3

doLoop1:
sra r6, 1 ; l >>= 1
mv  r14, r4
add r14, r6 ; r14 = mr+l
mv  r11, r14  ; ***********************************

mv  r14, r5
sub r14, r11 ; r14 = nn - (mr+l)
blz r14, doLoop1 ; while( mr+l > nn )

mv  r14, r6
addi r14, -1
and r14, r4
mv  r4, r14
add r4, r6 ; mr = (mr & (l-1)) + l

mv  r14, r4
sub r14, r2 ; r14 = mr - m
blez r14, contditForLoop ; if(mr <= m) continue;

mv  r7, r1
addi r7, 1 ; r7 = fr

mv  r8, r7
add r8, r3 ; r8 = fi

mv  r14, r7
add r14, r2 ; r14 = &fr[m]
load r13, 0(r14) ; r13 = fr[m]
mv  r12, r7
add r12, r4 ; r12 = &fr[mr]
load r11, 0(r12) ; r11 = fr[mr]
store r11, 0(r14) ; fr[m] = fr[mr]
store r13, 0(r12) ; fr[mr] = fr[m]
mv  r14, r8
add r14, r2 ; r14 = &fi[m]
load r13, 0(r14) ; r13 = ti = fi[m]
mv  r12, r8
add r12, r4 ; r12 = &fi[mr]
load r11, 0(r12) ; r11 = fi[mr]
store r11, 0(r14)
store r13, 0(r12)

contditForLoop:
addi r2, 1 ; ++m
bez r0, ditForLoop ; end ditForLoop
endDitForLoop:

; recycle r4, r5

li r6, 1 ; l = 1
load r4, 0(r1)
addi r4, -1 ; r4 = k = m - 1

whileMain:
mv  r14, r3
sub r14, r6
blez r14, endWhileMain ; while (l < n) {

mv  r5, r6
sll r5, 1 ; r5 = istep = l << 1;

li r2, 0 ; m = 0
outerForLoop: ; for (m=0; m < l; ++m)
mv  r14, r6
sub r14, r2
blez r14, endOuterForLoop


mv r9, r2 ; r9 = j = m

sllv r9, r4

;mv r14, r4; r14 = temp = k
; n2loop:
; bez r14, n2loopdone
; sll r9, 1
; addi r14, -1
; bez r0, n2loop;
; n2loopdone: ; r9 = j = m << k;

; load twiddle factors
mv  r12, r8
add r12, r3 ; r12 = sinewave = fi + n
mv  r14, r3
sra r14, 2 ; n/4

add r14, r9, ; j + n/4
add r14, r12 ; &sinewave[j+n/4]
load r10, 0(r14) ; r10 = wr = sinewave[j+n/4]
mv  r14, r12
add r14, r9 ; &sinewave[j]
load r11, 0(r14) ; r11 = wi = sinewave[j]

mv r13, r11
li r11, 0
sub r11, r13 ; wi = -wi
sra r10, 1 ; wr >>= 1
sra r11, 1 ; wi >>= 1

; free up some registers, forget r1
store r4, 3(r0) ; store k
store r5, 4(r0) ; store istep
store r3, 2(r0) ; store n
store r6, 5(r0) ; store l


mv r1, r2 ; r1 = i = m
store r2, 1(r0) ; store m

innerForLoop:
load r3, 2(r0) ; load n
mv  r15, r3
sub r15, r1 ; n - i

blez r15, endInnerForLoop ; for (i=m; i<n; i+=istep) {

load r6, 5(r0) ; load l
mv  r9, r1
add r9, r6 ; j = i + l

mv r12, r10 ; r12 = wr
mv  r15, r7
add r15, r9 ; &fr[j]
load r13, 0(r15) ; r13 = fr[j]
li r6, 1
jalr r6, r6 ; call FIX_MPY
mv r2, r12  ; tr = FIX_MPY(wr,fr[j])

mv r12, r11 ; r12 = wi
mv  r15, r8
add r15, r9 ; &fi[j]
load r13, 0(r15) ; r13 = fi[j]
li r6, 1
jalr r6, r6 ; call FIX_MPY
sub r2, r12 ; tr = FIX_MPY(wr,fr[j]) - FIX_MPY(wi,fi[j]);

mv r12, r10 ; r12 = wr
mv  r15, r8
add r15, r9 ; &fi[j]
load r13, 0(r15) ; r13 = fi[j]
li r6, 1
jalr r6, r6 ; call FIX_MPY
mv r3, r12 ; ti = FIX_MPY(wr,fi[j])

mv r12, r11  ; r12 = wi
mv  r15, r7
add r15, r9 ; &fr[j]
load r13, 0(r15) ; r13 = fr[j]
li r6, 1
jalr r6, r6 ; call FIX_MPY
add r3, r12 ; ti = FIX_MPY(wr,fi[j]) + FIX_MPY(wi,fi[j]);
mv  r15, r7
add r15, r1 ; &fr[i]
load r4, 0(r15) ; qr = fr[i]
mv  r15, r8
add r15, r1 ;
load r5, 0(r15) ; qi = fi[i]

sra r4, 1 ; qr >>= 1
sra r5, 1 ; qi >>= 1

mv  r14, r4
sub r14, r2
mv  r15, r7
add r15, r9
store r14, 0(r15) ; fr[j] = qr - tr

mv  r14, r5
sub r14, r3
mv  r15, r8
add r15, r9
store r14, 0(r15) ; fi[j] = qi - ti

mv  r14, r4
add r14, r2
mv  r15, r7
add r15, r1
store r14, 0(r15) ; fr[i] = qr + tr

mv  r14, r5
add r14, r3
mv  r15, r8
add r15, r1
store r14, 0(r15) ; fi[i] = qi + ti


load r5, 4(r0) ; load istep
add r1, r5 ; i += istep
bez r0, innerForLoop
endInnerForLoop:

load r4, 3(r0) ; load k
load r5, 4(r0) ; load istep
load r3, 2(r0) ; load n
load r6, 5(r0) ; load l
load r2, 1(r0) ; load m
addi r2, 1 ; m++
bez r0, outerForLoop
endOuterForLoop:

addi r4, -1 ; --k
mv r6, r5  ; l = istep
bez r0, whileMain

endWhileMain: